loading the R packages
## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 1.0.0
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
##
##
## Attaching package: 'plotly'
##
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
##
## The following object is masked from 'package:stats':
##
## filter
##
##
## The following object is masked from 'package:graphics':
##
## layout
##
##
## Please note that rgdal will be retired during 2023,
## plan transition to sf/stars/terra functions using GDAL and PROJ
## at your earliest convenience.
## See https://r-spatial.org/r/2022/04/12/evolution.html and https://github.com/r-spatial/evolution
## rgdal: version: 1.6-3, (SVN revision 1187)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.5.2, released 2022/09/02
## Path to GDAL shared files: C:/Users/hajam/AppData/Local/R/win-library/4.2/rgdal/gdal
## GDAL binary built with GEOS: TRUE
## Loaded PROJ runtime: Rel. 8.2.1, January 1st, 2022, [PJ_VERSION: 821]
## Path to PROJ shared files: C:/Users/hajam/AppData/Local/R/win-library/4.2/rgdal/proj
## PROJ CDN enabled: FALSE
## Linking to sp version:1.5-1
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading sp or rgdal.
##
## Loading required package: spData
##
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
##
## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired during 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
##
## rgeos version: 0.6-1, (SVN revision 692)
## GEOS runtime version: 3.9.3-CAPI-1.14.3
## Please note that rgeos will be retired during 2023,
## plan transition to sf functions using GEOS at your earliest convenience.
## GEOS using OverlayNG
## Linking to sp version: 1.5-1
## Polygon checking: TRUE
##
##
##
## Attaching package: 'transformr'
##
##
## The following object is masked from 'package:sf':
##
## st_normalize
##
##
## Loading required package: Matrix
##
##
## Attaching package: 'Matrix'
##
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
##
## Loading required package: foreach
##
##
## Attaching package: 'foreach'
##
##
## The following objects are masked from 'package:purrr':
##
## accumulate, when
##
##
## Loading required package: parallel
##
## This is INLA_23.01.02-1 built 2023-01-02 19:28:15 UTC.
## - See www.r-inla.org/contact-us for how to get help.
You can also embed plots, for example:
load("perfecture.rds")
glimpse(reduced_dataset)
## Rows: 1,128
## Columns: 5
## Groups: birth_year [24]
## $ birth_year <dbl> 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995, 1995,…
## $ prefecture <chr> "01", "02", "03", "04", "05", "06", "07", "08", "09", "10",…
## $ lbw <dbl> 4087, 998, 887, 1604, 718, 783, 1652, 2206, 1681, 1383, 520…
## $ normal <dbl> 45866, 12976, 12133, 20667, 9279, 10726, 19657, 26056, 1698…
## $ total <dbl> 49953, 13974, 13020, 22271, 9997, 11509, 21309, 28262, 1866…
map <- read_sf("jpn_adm/jpn_admbnda_adm1_2019.shp")
ggplot(map) + geom_sf()
map$ADM1_PCODE <- gsub( "JP", "", as.character(map$ADM1_PCODE))
map <- map %>% select(ADM1_EN, ADM1_PCODE, geometry)
data_japan <- reduced_dataset %>% rename(ADM1_PCODE = prefecture)
d <- data_japan %>% select(ADM1_PCODE, birth_year, lbw)
data_japan <- data_japan[order(data_japan$ADM1_PCODE, data_japan$birth_year),]
checking the first several row
data_japan <- data_japan %>% select(ADM1_PCODE, birth_year, lbw, total)
data_japan[1:20,]
## # A tibble: 20 × 4
## # Groups: birth_year [20]
## ADM1_PCODE birth_year lbw total
## <chr> <dbl> <dbl> <dbl>
## 1 01 1995 4087 49953
## 2 01 1996 4079 49794
## 3 01 1997 4108 48924
## 4 01 1998 4333 49065
## 5 01 1999 4217 46691
## 6 01 2000 4349 46791
## 7 01 2001 4224 46253
## 8 01 2002 4406 46122
## 9 01 2003 4169 44948
## 10 01 2004 4275 44029
## 11 01 2005 4073 41421
## 12 01 2006 4280 42206
## 13 01 2007 4150 41557
## 14 01 2008 3977 41077
## 15 01 2009 3726 40164
## 16 01 2010 3988 40168
## 17 01 2011 3847 39308
## 18 01 2012 3810 38692
## 19 01 2013 3780 38195
## 20 01 2014 3645 37068
library(SpatialEpi)
n.strata <- 1
E <- expected(population = data_japan$total, data_japan$lbw, n.strata = 1)
Creating a dataframe for expected counts
nyears <- length(unique(data_japan$birth_year))
admi1E <- rep(unique(data_japan$ADM1_PCODE), each = nyears)
ncountries <- length(unique(data_japan$ADM1_PCODE))
yearsE <- rep(unique(data_japan$birth_year), times = ncountries)
JapanE <- data.frame(ADM1_PCODE = admi1E, birth_year = yearsE, E = E)
head(JapanE)
## ADM1_PCODE birth_year E
## 1 01 1995 4584.754
## 2 01 1996 4570.161
## 3 01 1997 4490.311
## 4 01 1998 4503.252
## 5 01 1999 4285.363
## 6 01 2000 4294.541
d <- full_join(data_japan, JapanE, by = c("ADM1_PCODE", "birth_year"))
calculating SIR which is the observed divided by the expected
d$SIR <- d$lbw/d$E
SIR_data <- d %>% select (ADM1_PCODE, birth_year, SIR)
head(SIR_data)
## # A tibble: 6 × 3
## # Groups: birth_year [6]
## ADM1_PCODE birth_year SIR
## <chr> <dbl> <dbl>
## 1 01 1995 0.891
## 2 01 1996 0.893
## 3 01 1997 0.915
## 4 01 1998 0.962
## 5 01 1999 0.984
## 6 01 2000 1.01
names_of <- read_delim("prefectures.txt", delim = ",")
## Rows: 47 Columns: 4
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): id, prefecture_id, prefecture_en, prefecture_ja
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
names_of <- names_of %>% select(prefecture_id, prefecture_en) %>% rename(ADM1_PCODE = prefecture_id)
SIR_data <- full_join(SIR_data, names_of, by = "ADM1_PCODE")
save(SIR_data, file = "SIR_data.Rdata")
g <- ggplot(SIR_data, aes(x = birth_year, y = SIR,
group = prefecture_en, color = prefecture_en)) +
geom_line() + geom_point(size = 2) + theme_bw()
g
ggplotly(g)
map <- full_join(SIR_data, map, by = "ADM1_PCODE")
map <- st_as_sf(map)
ggplot(map) + geom_sf(aes(fill = SIR)) +
theme_bw() +
theme(
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank()
) +
scale_fill_gradient2(
midpoint = 1, low = "blue", mid = "white", high = "red"
) +
transition_time(birth_year) +
labs(title = "Year: {round(frame_time, 0)}")
load("japanstuartv2.rdata")
head(pnb)
## [[1]]
## [1] 0
##
## [[2]]
## [1] 3 5
##
## [[3]]
## [1] 2 4 5
##
## [[4]]
## [1] 3 5 6 7
##
## [[5]]
## [1] 2 3 4 6
##
## [[6]]
## [1] 4 5 7 15
nb2INLA("map.adj", pnb)
g <- inla.read.graph(filename = "map.adj")
d$idarea <- as.numeric(as.factor(d$ADM1_PCODE))
d$idarea1 <- d$idarea
d$idtime <- 1 + d$birth_year - min(d$birth_year)
formula <- d$lbw ~ f(idarea, model = "bym", graph = g) +
f(idarea1, idtime, model = "iid") + idtime
res <- inla(formula,
family = "poisson", data = d, E = E,
control.predictor = list(compute = TRUE)
)
d$RR <- res$summary.fitted.values[, "mean"]
d$LL <- res$summary.fitted.values[, "0.025quant"]
d$UL <- res$summary.fitted.values[, "0.975quant"]
d$logrr <- log(d$RR)
d$exprr <- exp(d$RR)
d$ADM1_PCODE <- as.numeric(as.character(d$ADM1_PCODE))
map$ADM1_PCODE <- as.numeric(as.character(map$ADM1_PCODE))
map <- full_join(map, d, by = c("ADM1_PCODE", "birth_year"))
ggplot(map) + geom_sf(aes(fill = RR)) +
theme_bw() +
theme(
axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank()
) +
scale_fill_gradient2(
midpoint = 1, low = "blue", mid = "white", high = "red"
) +
transition_time(birth_year) +
labs(title = "Year: {round(frame_time, 0)}")